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Abstract. The linear stability of isothermal Bondi accre- 
tion with a shock is studied analytically in the asymptotic 
limit of high incident Mach number M. i . The flow is un- 
stable with respect to radial perturbations as expected 
by Nakayama (1993), due to post-shock acceleration. Its 
growth time scales like the advection time from the shock 
r s h to the sonic point r son . The growth rate of non-radial 

2 /3 

perturbations / = 1 is higher by a factor M.^ , and is 
therefore intermediate between the advection and acoustic 
frequencies. Besides these instabilities based on post-shock 
acceleration, our study revealed another generic mecha- 
nism based on the cycle of acoustic and vortical pertur- 
bations between the shock and the sonic radius, indepen- 
dently of the sign of post-shock acceleration. The vortical- 
acoustic instability is fundamentally non-radial. It is fed 
by the efficient excitation of vorticity waves by the isother- 
mal shock perturbed by acoustic waves. The growth rate 
exceeds the advection rate by a factor log .Mi. Unstable 
modes cover a wide range of frequencies from the funda- 
mental acoustic frequency ~ c/r s h up to a cut-off ~ c/r son 
associated with the sonic radius. The highest growth rate 
is reached for I = 1 modes near the cut-off. The addi- 
tional cycle of acoustic waves between the shock and the 
sonic radius is responsible for variations of the growth 
rate by a factor up to 3 depending on its phase relative 
to the vortical-acoustic cycle. The instability also exists, 
with a similar growth rate, below the fundamental acous- 
tic frequency down to the advection frequency, as vorticity 
waves are efficiently coupled to the region of pseudosound. 
These results open new perspectives to address the stabil- 
ity of shocked accretion flows. 

Key words: Accretion, accretion disks - Hydrodynamics 
- Instabilities - Shock waves - Binaries: close - X-rays: 
stars 



1. Introduction 

Hydrodynamic instabilities in accretion flows can help un- 
derstanding the variability observed in the luminosity of 
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X-ray binaries. Numerical simulations have revealed the 
existence of such a hydrodynamical instability in the ac- 
cretion flow of a gas on a compact accretor moving at 
supersonic velocity (Bondi-Hoyle-Lyttleton accretion). A 
first step towards understanding the physical mechanism 
underlying this instability was made by Foglizzo & Tagger 
(2000, hereafter FTOO) who recognized the unstable cycle 
of entropic and acoustic perturbations between the shock 
and the sonic surface. This cycle is unstable if there is 
a large enough temperature difference between the shock 
and the sonic surface. The academic case of shocked Bondi 
accretion was studied by Foglizzo (2001, hereafter F01) 
who revealed the importance of non radial perturbations, 
and vorticity in particular. Both vorticity and entropy 
perturbations are advected from the shock to the accre- 
tor, and both are coupled to the acoustic perturbations. 
This coupling was formulated in a compact way by Howe 
(1975). If the adiabatic index is in the range 1 < 7 < 5/3, 
entropy and vorticity perturbations are intimately related 
in the shocked Bondi flow (Foglizzo 2002, in preparation), 
so that it is difficult to identify their respective roles in 
the instability mechanism. By contrast, their roles are 
well separated in the isothermal limit (7 = 1), where 
entropy perturbations are absent from the problem. The 
present paper is therefore dedicated to the study of the lin- 
ear stability of shocked accretion in the isothermal Bondi 
flow, where the incident Mach number is the only param- 
eter. The stability of shocked isothermal flows was stud- 
ied by Nakayama (1992, 1993) in the more general con- 
text of flows with small angular momentum. The study of 
Nakayama was restricted to axisymmetric perturbations, 
thus precluding any possible vortical acoustic cycle. In this 
approximation, Nakayama analytically obtained the result 
that the flow is unstable if the flow accelerates immedi- 
ately after the shock surface. In this respect the shocked 
Bondi flow should be unstable. We are therefore inter- 
ested here in an extension of Nakayama's results to the 
case of non radial perturbations, in order to include the 
effect of vorticity. It should be noted that isothermal flows 
are very particular as far as Bondi-Hoyle-Lyttleton (here- 
after BHL) accretion is concerned, with unsettled issues 
concerning the influence of the numerical resolution. The 
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Fig. 1. Radial profile of the Mach number in the unper- 
turbed shocked flow, for Mi = 5 (full line). The subsonic 
cavity stands between the sonic point and the shock ra- 
dius. The dashed lines show the analytical continuation of 
the solutions beyond the shock radius. 

instability observed in numerical simulations seems very 
weak in 3-D according to Ruffert (1996), whereas it is vio- 
lent in 2-D (Ishii et al. 1993, Shima et al. 1998). Pogorelov 
et al. (2000) discussed the possible responsability of the 
numerical procedure in producing the instability. With a 
different approach, the instability of the shocked Bondi 
flow described in this paper could contribute to guide our 
physical understanding of more complex flows involving 
isothermal shocks. 

The paper is organized as follows. Perturbed equations 
are described in Sect. || and eigenfrequencies are deter- 
mined numerically in Sect. ||. Analytical methods are used 
to disentangle the effect of advection from the effect of 
the boundary. In the spirit of F01, the coupling between 
the vorticity and acoustic perturbations in the classical 
isothermal Bondi flow, without a shock, is described in 
Sect. H Boundary conditions are taken into account to 
build a vortical acoustic-cycle in Sects. || and [| The insta- 
bility due to post-shock acceleration is analysed in Sect. [?]. 
The relationship between the vortical-acoustic instability 
and existing numerical simulations of BHL accretion and 
shocked discs is discussed in Sect. |8[ 

2. Linearized equations of the shocked Bondi flow 

2.1. Properties of the unperturbed shocked flow 

The Mach number in the unperturbed transonic flow satis- 
fies the following equation deduced from the conservation 
of the mass flux and Bernoulli constant B and the regu- 
larity at the sonic radius: 



The sonic radius is half of the Bondi radius GM/c 2 . The 
Bernoulli constant, 

B^ + cHogp-^, (2) 

is conserved along the flow lines. By contrast with the adi- 
abatic case 7 > 1, this quantity is not conserved through a 
shock. The shock radius r s h corresponding to an incident 
Mach number Mi is determined from Eq. (|l]) together 
with the Rankinc Hugoniot jump condition Mi = \j M.\. 
The Mach number profile is shown in Fig. |l|. It should be 
noted that the presence of a radial shock is not guaran- 
teed a priori, since the supersonic preshock flow (r > r sn ) 
could be continued without a shock down to the accre- 
tor (dashed line in Fig. [[]). Several physical mechanisms 
could trigger the formation of a shock, such as the heat- 
ing of protons to temperatures at which the fluid becomes 
collisionless (Mezaros & Ostriker 1983), the trapping of 
relativistic particules (Protheroe & Kazanas (1983) or the 
dissipation of magnetic fields (McCrea 1956, Scharlemann 
1981). In the context of Bondi-Hoyle-Lyttleton accretion, 
the existence of a stagnation point behind the accretor 
implies that a fraction of the supersonic flow decelerates 
to subsonic velocities, and therefore naturally produces a 
shock. The shock radius deduced from Eq. (|l|) increases 
with Mi for spherical accretion: 

^-elMf. (3) 

' son 

Compared to BHL accretion, the presence of an accre- 
tion shock at a distance exceeding the accretion radius 
2GM/v 2 XD is rather artificial and due simply to the as- 
sumption of purely radial velocity. In the highly super- 
sonic limit Afi > 1 initially considered by Hoyle & Lyt- 
tleton (1939), any trajectory with an angular momentum 
larger than ~ 2GM/voo would of course miss the accre- 
tor. Bearing this feature in mind, the case of strong shocks 
studied in this paper is still very useful in order to under- 
stand the mechanisms involved, because of the separation 
of timescales it enables. 

2.2. Scaling of timescales in the shocked Bondi flow 

If the shock is strong, the timescale r ac i v associated with 
advection is much longer than the sound crossing time 
r s h/c: 

dr Mi r sh 2r sh 
Tadv = / r-r ~ -= > ■ 4 

Jr Boa M 3 C C 

The third fundamental timescale of the flow is related to 
the presence of the sonic radius ~ r son /c. Together with 
Eq. (H), the scaling of the three special timescales of the 
problem is thus the following for strong shocks: 
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Fig. 2. Schematic structure of the perturbed acoustic (full 
line) and vorticity (dotted line) fields, depending on the 
frequency of the perturbation. The curve beyond the shock 
radius shows the acoustic perturbation in the Bondi flow 
without a shock. The turning point r t separates the region 
of pseudosound from the region of propagation of acoustic 
pertrbations. The vortical-acoustic coupling is most effi- 
cient in the region between r son and r c g . 

2. 3. Differential system for the perturbations of the 
Bondi flow 

As in the case 7 > 1 studied by F01, the vorticity equation 
can be directly integrated, so that the three quantities 
(r 2 5w r , rvSwg, rv5w v ) are conserved when advected. The 
perturbations of the radial velocity Sv r and density dp in 
the isothermal Bondi flow are conveniently described by 
the functions /, g defined by: 



2 S P 



f — vSv r + c — , 
P 

_ Sv r dp 



(6) 
(7) 



They satisfy the same differential equations as the func- 
tions f,g in the Appendix B of F01, by suppressing the 
entropy terms: 



iuiv 2 g 



df_ imM 2 f 

V ~dr + 1 -M 2 ~ 1-M 2 ' 

dg iu>A4 2 g i^>p 2 f 

V dr~ + 1 -M 2 = c 2 (l -M 2 ) ^ ~ 



i5K s h iu> f r 
e Jr s 



where the function /i(r, u, I) is defined by: 



1-^ 

2 ' 



(8) 
(9) 

(10) 
(11) 



The constant 8K is related to vorticity as follows 

5K = r 2 v • (V x w). (12) 



By contrast with the case of adiabatic flows (F01), SK 
is the unique source term for the excitation of acoustic 
waves. For radial perturbations (I = 0, SK = 0), the equa- 
tions analysed by Nakayama (1992, 1993), with "J = f /iuj, 
are recovered. 

The position of the turning point rt of non-radial acoustic 
waves is defined by /1 = in Eq. (|l^). This turning point 
lies inside the subsonic cavity if u> > w s h, with 



LU sh = I? (I + 1)- 



(1-M 2 h )t—K—. 

r s h r sh 



(13) 



Below o; s h (or inside the turning point rt if it exists) , the 
acoustic perturbation is not propagating. It was named 
"pseudosound" by Ffowcs Williams (1969). 
The threshold between absorbed and trapped sound was 
introduced in F01 through the cut-off frequency uf ut , 
which is independent of the shock strength: 



1 



UJ, 



(14) 



This cut-off frequency used to be refered to as a "refrac- 
tion" cut-off in F01 for 7 > 1, due to the non- homogeneity 
of the sound speed which bends trajectories outwards. By 
contrast, in the isothermal flow considered here, acoustic 
trajectories are only bent inwards, due to the sole effect of 
flow acceleration. The existence of a cut-off for non-radial 
acoustic waves can be understood by considering the fre- 
quency dependence of the direction ip of propagation of 
the incoming acoustic flux relative to the radial direction. 
This angle is approximated in Appendix E in the WKB 
approximation (Eq. E.ll). The higher the frequency, the 
smaller -0, so that high frequency waves with a given or- 
der / are pointing towards the accretor and are absorbed 
inside the sonic sphere. 

The scaling of frequencies stressed in Eq. (|J) for strong 
shocks consequently separates three ranges of eigenmodes 
schematized in Fig. || 

(i) absorbed sound ujf ut < u>, 

(ii) trapped sound w s h < to < W; cut , 

(iii) pseudosound l/r a dv < uj < lo^. The correspond- 
ing structure of the vorticity wave is also displayed in 
Fig. |^, in each range of frequencies. 

2.4- Boundary condition at the shock 

The boundary conditions at the shock associated with the 
variables /, g can be computed with the same method as 
Nakayama (1992), extended to the case of non radial per- 
turbations. Let us denote by AC the radial displacement of 
the shock produced by a sound wave propagating against 
the stream in the subsonic region , and Av — —iuA^ its 
velocity. The incident Mach number M.\ in the frame of 
the shock is altered by the perturbations of velocity and 
displacement as folllows: 

M[(r sh + AC) - A4i(r Bh ) + — + AC^i, (15) 

c or 
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-Mi(r sh ) 



1 



ir/cM^ 



Av 

1 

c 



(16) 



where the parameter 77 measures the strength of the local 
gradient of the M ach number immediately after the shock: 



n 



l-M? 



1 - 



' son 



(17) 



aiogr x-,v l2 

The boundary conditions / s h, <? s h, written as functions of 
Av, are deduced from the conservation of the mass flux 
and impulsion across the shock: 

Av 

2) ' ' " 



sh 



-tfM 2 (l -Mi) 
Av 



(1 



M 2 " 



Mi 



l'2 



(18) 



g sh = {I -M z 2 ) — . 

V2 



(19) 



The non radial perturbation of the velocity, and the per- 
turbed vorticity computed in Appendix A enable us to 
relate SK s h to Av: 



5K sh = l(l + l)c 2 (l 



Ml? 



1 



irjcMi 



Av 

V2 ' 



(20) 



Eq. (|2(]) can be used together with Eqs. (p^-^9|) in order 
to express / s h, <7sh as functions of 5K s h for non-radial per- 
turbations. In Appendix C the boundary value problem is 
reduced to a single equation incorporating the boundary 
conditions both at the shock and at the sonic point. The 
eigenfrequencies satisfy the following equation, where /o is 
the unique homogeneous solution which is regular at the 
sonic point: 



d 2 fo 
dr 2 



d log f 1 — M' 

dr 



M 



dfo 
dr 



2iuM 
(1-M 2 )c 

4^f /o = °' (21) 



uvr sh 



(l-M 2 2 )v-M. 
1 + Ml iujr sh 



iujr sh 



1(1 + 1) 



c 



1 - 
M2V 



M 2 



2 c 

r a h f 
JO 



M2V 



dfo , v 
fo(r sh ) _ 



r s h 

, r r 1+M 2 d, 



To, 



M 



dr. (22) 



Eq. (|22j) is independent of the normalization of fo- It 
should be noted that the integral on the right hand side 
is well defined despite the singularity of the phase near 
the sonic radius. This equation describes the perturbation 
of the shock by the interplay of the acoustic and vortical 
perturbations. The same calculation in any other poten- 
tial (e.g. the Paczynski-Wiita potential) would lead to the 
same system of equations (21 -p2|), only the shape of M 
(and its derivative rj) described by Eq. (Q) would be af- 
fected. The left hand side of Eq. ( p2| ) involves only the 
acoustic perturbation fo, and is independent of the vor- 
tical perturbation. By contrast, the integral on the right 
hand side describes the acoustic feed back of the vortical 
perturbation coupled to the acoustic field. 



— 1 = 

— Nakayama (1993) 

1 - 1 max 

I = 1 high 

1 = 1 low freq. 




Fig. 3. Growth rate of the purely growing radial instabil- 
ity, multiplied by the advection time (thick dashed line). 
The thin dashed lines correspond to the lower and upper 
bounds determined by Nakayama (1993). Anticipating on 
the calculations of the rest of the paper, the four other 
curves are analytical estimates of the growth rate of I = 1 
perturbations at high Mach number. The growth rate of 
the entropic-acoustic instability increases along the dot- 
ted line at low frequency (Eq. |74|) , and is bounded by the 
two thin lines at high frequency (Eqs. ^37|^). The thick 
full line describes the most unstable branch (Eq. |85|). 



3. Spectrum of eigenfrequencies 

3.1. Numerical procedure to determine the spectrum of 
eigenfrequencies 

The regularity of the solution at the sonic radius has al- 
ready been discussed in F01: the sonic point is a regular 
singularity if 7 < 5/3. For a given incident Mach number 
.Mi, numerical intergration is performed from the sonic 
point towards the shock in order to determine the eigen- 
frequencies of perturbations with a latitudinal number I. 
For a given value of u>, the unique solution which is regu- 
lar at the sonic point is expanded in a Frobenius series in 
order to start the numerical integration away from the sin- 
gularity. A Runge-Kutta algorithm is then used to simul- 
taneously integrate four functions, namely ,fp(r s h), 5o(?"sh), 
the integral on the right hand side of Eq. (g2|) and the in- 
tegral phase inside it. w is varied and shooting is repeated 
until Eq. (|2^ ) is satisfied. 

3.2. Purely growing radial instability 

The radial instability found numerically in Fig. ^ is consis- 
tent with the result of Nakayama (1992, 1993), who found 
that the local acceleration of the flow immediately after 
the shock is a source of non-oscillating radial instability. 
The growth rate is shown as a function of the incident 
Mach number, in units of the advection time r a d v - The 
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Fig. 4. Real and imaginary parts of the eigenfrequencies 
of the modes I = 0, 1, 2, 3 in the isothermal flow, for M.% = 
100 (cr a dv/47rr s h ~ 2.8). The real part is measured in units 
of the cut-off frequency wj ut (lower axis) and advection 
frequency (upper axis). The imaginary part is normalized 
to the advection time. Analytical estimates of the growth 
rates of the modes I = 0, 1, 2 at low freqency are indicated 
using big symbols (Eqs. (78), ^3| and |8^ in Sects. || and [?]). 

lower and upper bounds found analytically by Nakayama 
(1993, Eq. 27) are also displayed in Fig. |: 



< oji < 



(23) 



l + Mf- 1 + Mi 

where the parameter v describes the acceleration of the 
flow and is proportionnal to r\: 

1 /<9$ 2c 2 ' 
v 2 V dr r sh 



(Ml 1) 



(24) 
(25) 



Note that Eq. ( J24| ) is corrected for a factor 2 for spherical 
flows, as noticed by Nakayama (1993). The lower and up- 
per bounds of Eq. ( p3| ) cover a wide range of frequencies 
from the advection frequency ~ 2/(3r a dv) up to the acous- 
tic frequency ~ 2c/r s h- In the present case of isothermal 
radial accretion, the growth rate computed numerically 
converges precisely towards the value 2/r a dv 

3.3. Non-radial oscillatory instabilities 

The effect of non radial perturbations (I > 1) involves an 
integral on the right hand side of Eq. (|2^) , which reflects 
the cumulative excitation of acoustic waves by the vor- 
ticity perturbations advected from the shock to the sonic 
radius. Fig. ^ shows the eigenspectrum of the isothermal 
Bondi flow for Mi = 100, I — 0, 1, 2, 3. This spectrum is 
characterized as follows: 

(i) The most unstable mode is asymmetric (I = 1) and 
well separated from the other eigenmodes. The growth 
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Fig. 5. Real and imaginary parts of the eigenfrequencies 
of the modes I = 1 in the isothermal flow, for different 
Mach numbers. Frequencies are multplied by the advec- 
tion time. The strong I = 1 instability due to postshock 
acceleration corresponds to the thick line in both plots. In 
the upper plot, the growth rate of the fastest of the mode 
I = is also displayed (thick dashed line). The analytical 
approximation of the fastest growth rates tuf , w™ ax of the 
I = 1 instabilities due to the entropic-acoustic cycle and to 
post-shock acceleration are shown as the thin dashed line 
(Eqs. 68) and the thin dotted line Eq. (|8|). In the bottom 
plot, the cut-off frequency w™', the minimum acoustic fre- 
quency uj s h, and the frequency w™ ax of the most unstable 
mode (Eq. S3) are indicated as references. 



rate of this low frequency mode is significantly larger than 
the advection time, nearly four times faster than the un- 
stable radial mode. The analytical estimates of the growth 
rates at low frequency correspond to Eqs. ( [77| - [78| ), (^-|||) 
and (H-H) in Sects. § and @. 

(ii) Non radial unstable modes are very numerous and 
cover a wide range of frequencies, from the advection fre- 
quency (~ fshAsh), up to the cut-off frequency (oc c/r son ). 
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Fig. 6. Efficiency | Q a d v | of the excitation of acoustic waves 
propagating against the stream by the advection of a vor- 
tex in an isothermal flow, for I = 1,2,3. The frequency 
is in units of the cut-off frequency a;" 1 '. The dotted lines 
correspond to the asymptotic behaviour Q a d v °c w 
at low frequency determined analytically in Eq. (f 



(2i-l)/3 



(iii) A striking feature of the eigenspectrum in Fig. [| 
is the apparent oscillation of the imaginary part as a func- 
tion of the real part of the frequency. This behaviour is 
explained in Sect. ^| on the basis of the vortical- acoustic 
cycle. 

The variation of the growth rate with the incident Mach 
number is shown in Fig. |5|. The higher the Mach number, 
the more unstable the non radial modes, while the radial 
instability saturates. The growth rate of the most unsta- 
ble I = 1 mode increases with the incident Mach number 
much faster than in the other modes. The real part of its 
complex frequency is intermediate between the advection 
and acoustic frequencies. It is computed analytically in 
Sect, g 

The rest of the paper is dedicated to understanding this 
apparently complicated spectrum, and to questioning the 
role of post-shock acceleration. 



4. The sound of vorticity 

4-1- Frequency dependence of the vortical-acoustic 
coupling 

Let us assume that a vorticity wave is advected from infin- 
ity towards the accretor: what is the sound produced by 
the advection of this vorticity ? This question is solved in 
Appendix D by computing the coupling coefficient Q a d v 
measuring the outgoing acoustic flux associated with 
a vorticity perturbation SK: 



Ma 

r :2 



^IQadvH^I 2 , 



(26) 



The relationship between the acoustic flux F± and the 
acoustic perturbation f± is deduced from the WKB ap- 
proximation in Appendix B: 



F + 



c 2 M l/: 



(27) 



Thus the complex efficiency Q a dv is defined at a radius 
R 3> c/ui, in the region of propagation of acoustic waves, 
as follows: 



Qadv(i?) 



(_a_yf_ 

KM) SK 



1 



2lLUC 



; <= Jr t 1- 



Mr 2 



+ dr 



(28) 
dr. (29) 



Q a dv is approximately independent of R in the WKB re- 
gion of propagation of acoustic waves, i.e. away from their 
turning point. This coefficient is approximated in Ap- 
pendix D in the low frequency limit c/R <C ui <C c/r son , 
i.e. below the cut-off frequency, using the function TC(l) 
defined from an integral of the Spherical Bessel function 
ji (Eq. [Dj: 



\Q 



l + l 21^1 

3 3 e t 



adv 



H{1) 



l-3...(2/ + l) 2 



^ uir son y 



(30) 



with H(l) ~ 0.45, H(2) - 0.33 and H(S) ~ 0.28. The 
direct integration of Eq. (B3) is compared to the approx- 
imation obtained in Eq. ( p0[ ) for I = 1, 2, 3 in Fig. ||. The 
asymptotic behaviour |Q a dv| ^ uj 21 ^ 1 ^ 3 is correctly re- 
produced. The scaling factor of this power law obtained 
analytically is too low by a factor 2 for I = 1 perturba- 
tions, which can be attributed to the roughness of the 
approximation. This calculation clearly indicates that the 
vortical-acoustic coupling is most efficient near the cut-off 
frequency. The same conclusion was reached in F01 con- 
cerning the efficiency of the entropic-acoustic coupling. 

4 ■ 2. Region of coupling 

This approximation enables us to determine the region 
where the coupling between the vorticity and the acoustic 
perturbations is most efficient. It occurs mainly in the 
subsonic region r son < r < r c g where the wavelength of 
the vorticity perturbation is largest (see Fig. ||): 



LO I — ~ 2tt. 



(31) 



A similar result was obtained in F01 concerning entropy 
perturbations for 7 = 5/3. Due to the differences in veloc- 
ity profiles in Eq. ([3l|), r c s oc ui~ 2 l 3 for 7 = 5/3 at high 
frequency (Eqs. E.6, E.ll and E.12 of F01), whereas 



r c s oc 



^ ujr son y 



(32) 



in the isothermal flow at low frequency. It is therefore 
important to note that the most efficient coupling at low 
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frequency lies in the region of pseudosound well within the 
turning point (r t ~ c/lu) of acoustic waves: 



^son < r c g < r t . 



(33) 



According to Eq. (31 ), the upper bound r c s of the coupling 
region for perturbations near the advection frequency co- 
incides with the shock radius. 



5. Vortical-acoustic instability at high frequency 

5.1. The vortical-acoustic cycle for o; s h <C u> <C ajf nt 

The formalism developed in FTOO for the cntropic- 
acoustic cycle within the WKB approximation can be 
transposed to the case of the vortical-acoustic cycle. This 
formalism requires the frequency to be high enough so 
that acoustic waves can be identified and their propaga- 
tion time measured. An initial perturbation of vorticity 
6K(to,r s h) at the shock, advected towards the accretor, 
triggers the excitation of an acoustic flux propagating 
outwards. As it reaches the shock surface, this acoustic 
flux induces new vorticity perturbations 8K(tr) + tq, r s h), 
where tq is the duration of this cycle. The global efficiency 
Q of the cycle is naturally 



Q 



SK(t + T Q ,r sh ) 
SK(t ,r sh ) 



(34) 



Q and tq depend a priori on the frequency oj r and spatial 
structure I of the perturbation considered. The growth (or 
damping) rate of the cycle is identified with the imaginary 
part of the eigcnfrcquency, given by: 



Wj = — log | Q|. 



(35) 



The timescale tq is dominated by the advection timescale 
-r adv if the shock is strong (M 2 < 1, r sh » r t ): 



1 dr 



7"adv ■ 



(36) 



The global efficiency Q can be decomposed in terms of 
the efficiencies of the coupling between F_ and SK during 
advection (Q a dv) and at the shock (Q s h), defined imme- 
diately after the shock according to Eq. (|28|) and: 



Qsh 



M 2 y SK sh 
M J /- ' 

Q = SadvQshC-^. 



(37) 
(38) 



The modulus of Q s h is directly interpreted as a coupling 
efficiency between F_ and jd/f | 2 : 



\SK sh \ 



Mo 



\Q sh \ 2 F. 



(39) 



The efficiency Q s h is computed in Appendix E, following 
a method introduced by D'lakov (1958) and Kontorovich 



(1958, 1959). Q s h can be approximated by a WKB analy- 
sis when the wavelength of the perturbation is small com- 
pared to the lengtliscale of the local gradients of the flow, 
by keeping only the first order terms in M 2 and c/ur, at 
high frequency and for strong shocks: 



21(1 + 1) 

mI 



1 - M 2 



rjc\ 



(40) 



This calculation proves that the flow reacceleration does 
not affect the vortical-acoustic coupling at the shock for 



> M 



(41) 



and that Q s h is approximately independent of frequency in 
the range w s h <C u> < cuf ut . As a consequence of Eqs. (37) 
and (|l0|), the vorticity is produced at the shock with an 
efficiency proportional to the shock strength: 

Sp- 
P 



SK sh tx Mif- oc Mi- 



(42) 



The global efficiency Q should consequently peak near the 
cut-off frequency, and drop abruptly above it. The fre- 
quency dependence of Q appears in Fig. @ (top picture, 
full line) , as deduced from Fig. || and Eq. ( ftcf ) in the WKB 
approximation. Thus the instability should be strongest 
near the cut-off frequency with the following scaling, ob- 
tained from Eqs. (O), (|(| and <M) for uj sh <C u <C w z cut : 



Id 



Mi 
Mo 



(43) 



where the scaling factor Mqi depends on I through both 
Qadv and Q sh : 



Mr 



4 

3e3 



[1 • 3 . . . (2Z + l)] 2 1 



[l(l + l)Y 



m 2 {i) 



(44) 



Moi oc diverges for I 3> 10, thus favouring the in- 
stability of low degree modes. According to Fig. ^, the 
extrapolation of Eq. (E3j) to ui ~ ojf ut gives a very rough 
upper bound of |Q|(w; ), particularly overestimlating it 
for high degree perturbations. Although Moi ~ 1.96 ex- 
ceeds M02 ~ 0.95 and M03 ~ 0.69, the actual efficiency 
|Q| is maximal for the mode I = 1, as checked in Fi g. p| 
by multiplying the curve |Q adv | (Fig. |) by |Q sh | (Eq. 
Using Eqs. (^) and (|36|), the growth rate oji at high fre- 
quency is deduced from Eq. ([43|): 



oj sh < oj r < w z cut , 



1 



1 



log 



2t : 



adv 



log 



2^ 



3 Mi 

2(21-1) 2(2-1) 

Mi 3 



Mr 



(45) 
(46) 

(47) 



Eq. (|47]) indicates that the range of unstable frequencies 
below the cut-off frequency gets narrower for I > 2. This is 
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also confirmed by the eigenspectrum obtained numerically 
for I = 1,2,3 in Fig. |]. Note that the efficiency |Q| was 
estimated using | Q a d v | which was computed for perturba- 
tions with a purely real frequency. Thus the estimate of the 
growth rate in Eq. (^) implicitly assumes that the imag- 
inary part is small compared to the real part. This is true 
since log .M 1 /T a d v <C ^/ cut - This calculation proves that 
the vortical-acoustic cycle is unstable for strong shocks, 
with a growth rate exceeding the growth rate of the ra- 
dial mode, and a mechanism which is independent of the 
sign of the local flow acceleration immediately after the 
shock. A more refined description of the entropic-acoustic 
cycle is developped in the next section in order to under- 
stand the apparent oscillations in the spectrum obtained 
numerically in Fig. || and Fig. 0. 

5.2. Additonnal contribution of the acoustic cycle 

5.2.1. Dispersion relation for the double cycle 

The role of the purely acoustic cycle had been anticipated 
in FTOO. It is characterized by a time scale t-r and a global 
efficiency 1Z. Let us show how the simultaneous existence 
of the two cycles can explain the apparent oscillation in 
the eigenspectrum of Fig. [|. The following extension of the 
analysis of FTOO is valid for both entropic-acoustic and 
vortical-acoustic cycles, by replacing "vortical" by "en- 
tropic" . 

A perturbation / is influenced by the two cycles as follows: 



(48) 



f(t) = Qf(t-T Q )+llf(t-T R ). 



A solution of the form f(t) oc exp(— itot) satisfies Eq. 
if the complex eigenfrequency (u> ri u>i) is a solution of the 
following dispersion equation (eq. 25 of FTOO): 



(49) 



This dispersion relation is recovered in Appendix F as 
a WKB approximation of the exact dispersion relation 
(p2|). The analysis of the dispersion relation (|^) in Ap- 
pendix G enables us to extract physical information from 
the complicated eigenspectrum in Fig. ^, such as the ratio 
of timescales tq/t-r, and the dimensionless efficiencies \ Q\ 
and \TZ\. 

5.2.2. Ratio of the two cycles timescales tq/t-r 

The timescale t-r of the acoustic cycle for strong shocks 
can be written approximately: 



rr s h 



dr 



TK ~J rt T^Mi- 

According to Eqs. ( |36| 
strong shock is simply: 



■ 2^1 oc Mf 
c 



(50) 



t-r 



Mi 
6 



and (|50|), the ratio tq/t-r for a 

(51) 



20- + 



+0 _ ♦ _ + _ ♦ o* o 



= 1 

= 2 
= 3 



(fl r / 0)'" 

Fig. 7. Number of eigenmodes in each oscillation of the 
eigenspectrum in Fig. f|, identified as the ratio of the 
timescales of the two cycles. For strong shocks, tq/t-r ~ 
Mi/6 

As shown in Appendix G, this ratio is the average number 
of eigenmodes in each oscillation of the eigenspectrum. As 
an illustration, this number is reported on Fig. ^ for the 
modes I — 1, 2, 3 in the case Mi = 100 corresponding to 
Fig. ^, in good agreement with Eq. (|5l"|). 

5.2.3. Efficiencies \Q\ and \K\ 

As noticed in FTOO, the acoustic cycle can be neglected 
(u>i ~ log|Q|/Tg) if the efficiencies and timescales of the 
two cycles are such that \a\ <C 1, with 



K 



(52) 



If a is not negligible, the acoustic cycle can contribute to 
either stabilize or destabilize the vortical-acoustic cycle. 
The most stabilizing effect of the acoustic cycle is the ef- 
fective reduction of Q by a factor 2 (see Appendix G). By 
contrast, its destabilizing contribution can be much larger 
comparatively, if the acoustic time is short compared to 
the advection time. This is the case for strong shocks. 
Depending on the relative phases of the two cycles, the 
growth rate can then cover the following range: 



i , \Q 

T Q 1 + I 



U < Ui < — log -!5L 

a | t q 1 - \o 



(53) 



Conversely, the values of the dimensionless physical pa- 
rameters |<2|,|7\L| can be extracted from the eigenspec- 
trum, as in Fig. [|, by measuring 

(i) the extremal values lu^,uj^ , 

(ii) the period Au> r of the oscillations 

(iii) the average number n of modes per period. 

As demonstrated in Appendix G, the timescale of the 
acoustic cycle is simply t-r = 27r/Aw r . Thus the oscil- 
lations are well resolved in the eigenspectrum only if the 
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+ -■- R„ 
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Fig. 8. Global efficiencies \Q\, \1Z\ asociated to non radial 
perturbations I = 1,2,3 deduced from the eigenspectrum 
in Fig. ^ in the framework of the double cycle model (dot- 
ted line), compared with the efficiencies deduced from the 
calculation of Q adv , Q s h and TZ a dv,TZ s h 

advection time is significantly longer than the acoustic 
time (i.e. Mi > 6). 

The values of \Q\, \7Z\ are determined using the following 
relations: 



cosh tt^- Q i 
1 21 = — 7T J ~At7< exp2mr- 



Au 
Au T 



\K\ = 

where 
Qi = 



cosh(n — l)n 

sinh tiit^—^ 
cosh(n - lj^ 



exp 27r 



Au 



Auj r 



(54) 
(55) 



(56) 
(57) 



As an illustration, Fig. |8| is obtained from the eigenspec- 
trum of Fig. H using Eqs. (|54| ) and (55). The accuracy 
of the measurements of Aw r , Aujj and Uj is of course of 



the order 1/n. As a check of consistency, the value of |Q| 
obtained from the product of |Q a dv| (Fig. ||) and \Q s h\ 



(Eq. 40) is also displayed, showing an excellent agreement 
except for the lowest frequency modes where WKB ap- 
proximation breaks down. It is remarkable that \ Q\ seems 
to be maximized at low frequency for I = 1 whereas it is 
maximized near the cut-off frequency u ; cut for I > 2. The 
efficiency \1Z\ of the acoustic cycle is bounded by one and 
decreases to zero above the cut-off frequency, as could be 
expected from FT00 and F01. 



rate 



5.2.4. WKB approximation of 1Z and growth 

The global efficiency 1Z can be decomposed in terms of 
the efficiencies IZ^dv and 7£ s h of the coupling between the 
acoustic fluxes F±, such that the global efficiency 1Z is 



TZ — IZadvTZshC 



(58) 



The efficiency 7Z a dv measures the outgoing acoustic flux 
F_ produced by the deviation of an ingoing acoustic flux 



adv 



f_ 

/+' 

im f sh 2 » „ dr 



F- = \K. 



adv 



(59) 

(60) 
(61) 



IZ&d-v is simply the limit when 7 — > 1 of the efficiency 
computed in F01: it is close to unity below the cut-off fre- 
quency cof ut , and decreases exponentially above this cut- 
off. Conversely, an acoustic flux F_ reaching a shock pro- 
duces a reflected ingoing acoustic flux F + with the effi- 
ciency 7£ sh : 

f+ 

S- 



TZsh = -r- , 



(62) 



\TZ sh \ 2 F„ 



(63) 



Keeping only the first order terms in M2 and c/wr, the 
complex efficiency is computed in Appendix E: 



sh 



-1 + 2M2 



(64) 



Thus the global efficiency \R\ should be close to unity be- 
low the cut-off frequency. The product |72. s h7^ a dv| is dis- 
played in Fig. U (bottom picture), in good agreement with 
the value of \1Z\ deduced from Fig. ^ and Eq. (p5|). 
The extremal values lof of the growth rate of the vortical- 
acoustic cycle are computed in Appendix G (Eq. p.7| ) for 
strong shocks: 



1 , Igj 
— log 1 - | | 

TQ 1 =F \OL\ 



(65) 



The asymptotic scaling of a is deduced from Eqs. (ff3j), 
© and M) : 



a 



Mi 



log 



Mi 



Mqi V^ c 



\(2l-l) 



(66) 
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Together with Eq. (pq), and despite the roughness of the 
approximation in Eq. ( |30| ) near the cut-off frequency, the 
leading order of the asymptotic values of u>f is comparable 
to: 



1 



Ul; 



2T a dv 

3 



log 



Mi 

4:M l 



log- 



Mi 



2r a dv (9M ;)ilog3^ 



(67) 
(68) 



The values of Q, 1Z, t-r/tq are thus responsible for a dis- 
persion of the growth rate by a factor 3 in Eqs. ( |67| - 
|68| ) near the cut-off frequency, depending on the relative 
phases of the two cycles. This factor 3 is consistent with 
both Figs. | and § 



6. Vortical-acoustic instability at low frequency 

The region of efficient coupling between the vorticity and 
acoustic perturbations lies in the region of pseudosound 



according to Eq. (33). Thus it seems natural to expect an 



instability at low frequency as long as the coupling region 



determined from Eq. (32) lies inside the shock radius, i.e. 



down to the advection frequency. This case is similar to 
the hole tone instability (e.g. whistling kettle) studied by 
Chanaud & Powell (1965), or in the oscillations of im- 
pinging shear layers (see a review by Rockwell 1983). In 
some of these cycles, the vorticity perturbations are cou- 
pled to the acoustic field in the region of pseudosound. 
As stressed by Chanaud & Powell (1965), this does not 
preclude the use of the term "acoustic feedback" so that 
we may talk of a vortical-acoustic cycle even in this range 
of frequencies. Although the problem cannot be treated in 
the WKB limit, the homogeneous solution can be approx- 
imated by a Spherical Bessel function of the first kind in 
the region of pseudosound far from the sonic point. The 
calculation in Appendix H shows that the case I = 1 and 
I > 2 must be treated separately. In the domain of pseu- 
dosound 1 <C |^T a dv| "C M%, Eq. (|||) is reduced to: 



ILOT 



y 

M2C dr 2 



l(l + l)(l + r)-2)f 
'l + A 



9le 1 ' 



T 



(jWT a dv) 3 /o-(69) 



6.1. Low frequency global cycle I = 1 

The / = 1 acoustic perturbation is approximated in Ap- 
pendix H as fo oc ji(ujr/c), and Eq. ( |69| ) is transformed 
into 

^M\ (lLUT adv f = ^ + (iwrrfv)* T e—v (70) 

A branch of solutions with uj r S> u>i corresponds to a 
vortical-acoustic cycle between the shock and the sonic 



point: 



WiT adv - log 



+ - log |w r T adv |, 

if |tJT adv | <^m\, (71) 



log 



2 log All - - log |w r Tadv|, 

if |u;r adv | «Aif- (72) 



The maximum growth rate of this branch of solution is 
reached for a real frequency similar to the frequency of 
the local instability. 



imMi 

8 Mi 
9r adv ° g AV 



M = 



2f3§ 

5 *rl(|) 



1.17. 



(73) 
(74) 

(75) 



This growth rate is compared to the results of numerical 
calculations in Fig. [| This growth rate exceeds the growth 
rate of the vortical-acoustic instability at high frequency 
(Eq. ^) , but is asymptotically smaller than the maximum 
growth rate ujf reached when the purely acoustic and vor- 
tical acoustic cycles are in phase (Eq. |68[ ). 

6.2. Low frequency global cycle I > 2 

If I > 2, the acoustic feedback due to an isolated vortical 
perturbation near the shock is negligible compared to the 
integral effect for strong shocks. Eq. (p9) is reduced to: 



'l + A\ 2-1 
3r I I (iLUT &dv )— e— ^ =1-1. 



(76) 



The most unstable modes are therefore the modes / = 2, 
with 

2n7r 



uj r ~ , with 1 n < Mi 

7"adv 



log 3 

^adv 



> 0. 



(77) 
(78) 



These approximations are in excellent agreement with the 
numerical result of Fig. ^. 

By contrast with the case of high frequency perturbations, 
it seems difficult to separate completely the effect of post- 
shock acceleration from the global vortical-acoustic cycle 
in this range of frequencies. 

7. Effects of post-shock acceleration: a global 
instability with a local criterion 

can be rewritten for strong 



The boundary condition 
shocks as follows: 

fsh AC . . M c |r?| c 
— - 77 for M < — r . 

C r sh M2 r s h r son 



(79) 
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This equation is similar to the argument of Nobuta & 
Hanawa (1994) concerning the balance total pressures 
(thermal and dynamical) on both sides of the shock. 
Rather than treating the shock surface as a material sur- 
face pushed by the local pressures on both sides of it, 
Eq. (fn^) is interpreted as follows: an excess of Bernoulli 
perturbation / (which we may call " energy" ) on the sub- 
sonic side of the shock is associated with a displacement 
of the shock in the direction of increase of the local Mach 
number. If r/ < (Eq. [Tt]), an excess of / produces an 
outward displacement of the shock. This statement alone 
is not conclusive: instability occurs only if the flow is not 
able to evacuate this excess of energy. Although the dis- 
placement of the shock indeed liberates some potential 
energy locally, the instability depends on the leakage of 
energy through the sonic radius (Nakayama 1993). 

7.1. Asymptotic growth rate of the radial instability 

For radial perturbations (2 = 0), the only way of evacuat- 
ing the excess of energy is through acoustic perturbations. 
Acoustic energy, however, is trapped in the subsonic re- 
gion if the real frequency of the mode is low enough. Thus 
a low frequency instability is expected. Although the cri- 
terion of the instability is indeed local, its growth rate 
depends on the spatial structure of the acoustic perturba- 
tion from the sonic radius to the shock, as illustrated by 
Eq. (E2) for strong shocks: 



dlog/o 
<91ogr 



(r sh ) 



1 



(80) 



fo is approximated as the acoustic perturbation of a uni- 
form medium in spherical coordinates (see Appendix B), 
using a Spherical Bessel function fo oc jo(u)r/c) for r s h S> 
r son . The solution of the dispersion relation ( p0[ ) at low 
frequency is a purely imaginary eigenfrequency: 



o, 

-rj 3M 2 c 
3 + V r sh 
,d\v\ << 



-6- 



di 



Tadv 



(81) 
(82) 

(83) 



The fact that the growth time scales like the advection 
time is not obvious a priori. 



T) < 




n > o 




Fig. 9. Interaction of the vortex with the shock. The con- 
tribution of a vorticity perturbation 5K s h > to the 
Bernoulli perturbation is positive (vSv r > 0) , and induces 
an increase of the local accretion rate (g > 0). depend- 
ing on the sign of the postshock acceleration, a vorticity 
perturbation contributes to a local expansion (77 > 0) or 
collapse (77 < 0) of the shock. 

deduced from Eq. ( |l9| ) , we conclude that the vorticity con- 
tributes to the instability if v < 0. A quantitative calcu- 
lation of the growth rate requires taking into account the 
acoustic perturbation in the region from the shock to the 
sonic point. The fastest instability of the flow, well isolated 
from all the other modes of instability, is obtained by ne- 
glecting the last term on the right hand side of Eq. (|70|): 



5^32 



1 



2a r sh M % 

55 c 1 
2% r sh j^k 



1 ) 

3 



(84) 
(85) 



It is in excellent agreement with the numerical calcula- 
tions in Figs. |] and ||. This growth rate is intermediate 
between the advection and the fundamental acoustic fre- 
quencies. The real and imaginary parts being comparable, 

the growth rate is the time taken by the vorticity pertur- 

2 

bation to travel by one wavelength 27n;/w™ ax oc .Mf r s h <C 



7.2. Asymptotic growth rate of the I = 1 instability 

Non radial perturbations generate vorticity, which con- 
tribute to the energy balance in the subsonic part of the 
flow. Fig. H illustrates the contribution of a vorticity per- 
turbation to the Bernoulli perturbation through the ra- 
dial component of velocity v5v r . Thus the vorticity per- 
turbation participates to increase the local accretion rate 
(i.e. g > 0) in the regions where the shock moves inward 
(A£ < 0). Comparing this statement with g ~ iojA(/A4 2 c 



8. Discussion 

8.1. Relationship with other shock instabilities 

Most shock instabilities identified and studied in astro- 
physics are related to the acceleration or deceleration 
of the shock itself. The most famous is of course the 
Rayleigh- Taylor instability of accelerated shocks, anal- 
ysed by Bernstein & Book (1978). Decelerated shocks are 
also unstable to a rippling instability, studied by Vishniac 
(1983), Bertshinger (1986) and Vishniac & Ryu (1989). 
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By definition, stationary shocks are stable with respect 
to these mechanisms, but can still be unstable. Nakayama 
(1992, 1993) pointed out the radial instability of the shock 
if the flow is immediately accelerated after the shock, in 
isothermal flows. The validity of the postshock accelera- 
tion criterion for adiabatic flows is still uncertain, even 
for radial perturbations (Nakayama 1994). The vortical- 
acoustic cycle studied in the present paper resembles by 
many aspects to the entropic-acoustic cycle studied by 
FTOO and F01 in adiabatic flows. Two distinct mecha- 
nisms, however, are at work: 

(i) the vortical-acoustic cycle is fed by the vorticity 
production by a perturbed shock, which is highest for 
isothermal flows with a strong shock, 

(ii) the entropic-acoustic cycle is fed by the tempera- 
ture increase from the shock to the sonic radius, which is 
highest if 7= 5/3 (FTOO, F01). 

A calculation similar to Appendix E for adiabatic flows 
would show that the vorticity production by a perturbed 
shock is large only in the isothermal limit: 



SK 



sli 



/- 



1 

~M~2 



< 



27 
7-1 



(86) 



Although shocked adiabatic flows with 1 < 7 < 5/3 are 
subject to both entropic-acoustic and vortical-acoustic cy- 
cles in principle, their stability cannot be determined with- 
out a specific calculation. 



it is detached in high resolution simulations of adiabatic 
flows with 7 > 4/3. This can be understood qualitatively 
in terms of the weakness of pressure forces in isothermal 
flows compared to adiabatic flows (P oc p 1 ). More quanti- 
tatively, Foglizzo & Ruffert (1997) proved that the shock 
in the isothermal BHL flow cannot be detached unless 
the sonic surface extends up to distances comparable to 
the Bondi radius GM/c 2 , i.e. much larger than the accre- 
tion radius. In view of this topological difference between 
isothermal and adiabatic BHL flows, it seems easier to ex- 
tend the results obtained for shocked Bondi accretion to 
the subsonic region ahead of the accretor, for detached 
shocks, than in the very non-radial region of accretion 
behind the accretor for attached shocks. Although true, 
however, this argument is not conclusive. Even in a non 
radial isothermal flow, a small pressure perturbation of 
the shock is able to generate vorticity perturbations very 
efficiently. The lack of instability of isothermal BHL ac- 
cretion in 3-D must be sought for in the lack of acoustic 
feedback from the advected vorticity perturbation. Unfor- 
tunately we do not have quantitative arguments to explain 
why the acoustic feedback is so weak, apart from noticing 
that geometric compression of the vorticity perturbation 
from the shock to the sonic radius might be insufficient in 
the BHL flow. Any coherent explanation should of course 
also account for the strong instability observed in 2-D sim- 
ulations of BHL accretion (Shima et al. 1998). 



8. 2. Constraints on the instability mechanism of BHL 
accretion 

The stability of BHL accretion can be addressed with new 
tools, using the present results and the entropic-acoustic 
cycle of FTOO and F01. Although a detailed analysis of 
the existing numerical simulations of BHL accretion is 
postponed to a forthcoming paper, let us outline the pos- 
sible consequences of the entropic-acoustic and vortical- 
acoustic mechanisms on this specific flow. Numerical sim- 
ulations of BHL accretion in 3-D show a strong instability 
for adiabatic flows with 7 = 5/3 and small accretors (Ruf- 
fert & Arnett 1994, Ruffert 1994): this coincides nicely 
with the properties of the entropic-acoustic cycle, which is 
most unstable if the temperature gradient in the subsonic 
part of the flow is strongest (FTOO, F01). By contrast, 
nearly isothermal flows are rather stable in 3-D numeri- 
cal simulations (Ruffert 1996). This could seem puzzling 
given the strong vortical-acoustic instability described in 
the present paper: why would the entropic-acoustic cy- 
cle be relevant for the BHL instability with 7 = 5/3, 
and the vortical-acoustic cycle be irrelevant to the BHL 
stability for 7 = 1? This apparent contradiction can be 
partly solved by remembering the difference of topology 
between the spherical Bondi flow and the BHL flow, fol- 
lowing the remark of Sect. 



2.1. All the numerical simula- 



tions of isothermal BHL accretion seem to agree with the 
fact that the shock is attached to the accretor, whereas 



8.3. Relationship with the numerical simulations of 

shocked inviscid flows with low angular momentum 

Inviscid accretion flows with low angular momentum are 
much more complex than the Bondi flow in the sense that 
the shock position is not unique (Fukue 1987, Chakrabarti 
1989). The numerical simulations of Nobuta & Hanawa 
(1994) did not detect the effects of vorticity perturbations 
since they were restricted to axisymmetric motion. These 
simulations showed the evolution of the radial instabil- 
ity due to postshock acceleration: the unstable shock is 
either absorbed by the accretor or inflates to reach the 
stable outer position. 

Molteni, Toth & Kuznetsov (1999) performed non axisym- 
metric numerical simulations of shocked inviscid accretion 
flows with low angular momentum, for 7 = 4/3. Unstable 
m = 1 oscillations were observed, with no physical expla- 
nation. The entropic-acoustic and vortical-acoustic cycles 
might provide a physical basis to understand this insta- 
bility. 

8.4- Non linear evolution of the instability 

If the existence of the shock is postulated a priori, as in the 
present study, the flow might simply converge towards an- 
other solution, as illustrated by the dashed lines in Fig. [I]. 
The shock can either be absorbed by the accretor, leading 
to the fully supersonic solution, or could expand towards 
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infinity to establish the Bondi solution if the outer bound- 
ary condition allows it. 

In accretion flows stable with respect to the postshock 
acceleration criterion, but unstable through the vortical- 
acoustic (or entropic-acoustic) cycle, the instability might 
be saturated by the effect of the geometric dilution of 
the acoustic energy in the subsonic cavity. Let us assume 
that advected vorticity perturbations generate an acoustic 
flux propagating outward F_ <x |Q a dv| 2 |<^sh| 2 (Eq. p6| ), 
and that the acoustic perturbation produces in turn a 
vorticity perturbation l&Kstl oc M.i\6p-/p\ (Eq. f42|) . 
The amplitude of the pressure perturbation \5p—/p\ de- 
pends on the volume in which the acoustic flux F— is 

diluted through Eq. @: \8p-/p\ cx {r sh /r)F^ 2 . Thus 
the vortical-acoustic cycle is naturally stabilized when the 
shock reaches a distance r max defined by 



\Q\r s h, 



(87) 



where \Q\ is the global efficiency of the vortical- acoustic 
cycle in the linear regime. If the shock is not simply ab- 
sorbed by the accretor, the non linear evolution of the 
instability could lead to quasi periodic oscillations of am- 
plitude comparable to r max . However, given the number 
of unstable modes in a spectrum like Fig. ^, the vortical- 
acoustic instability might as well saturate into turbulence 
rather than be dominated by a single QPO. This issue can 
only be solved with numerical simulations. 

9. Conclusions 

The linear stability of shocked isothermal Bondi accre- 
tion has been studied by comparing the complex eigen- 
frequencies obtained through a direct numerical integra- 
tion (Sect. H]) to the analytical results obtained for strong 
shocks by two different methods: 

(i) an analytical estimate of the growth rate corre- 
sponding to a cycle of perturbations with a purely real 
frequency, obtained by separating the effects of advection 
from the boundary effects of the shock. This WKB ap- 
proximation, valid in the range of acoustic waves (c/r s h <C 
uj <C c/r son ), was used in Sects. |] and ||, 

(ii) an analytical estimate of the complex eigenfrequen- 
cies in the range of pseudosound (vshAsh *C u> <C c/r s h) 
using Spherical Bessel functions (Sects. || and |^). 

The results obtained by these methods are summarized as 
follows: 

- As expected by the postshock acceleration of 
Nakayama (1993), the isothermal Bondi accretion with a 
shock is unstable with respect to radial perturbations. Its 
growth rate is comparable to the advection time from the 
shock to the sonic point. 

- The present analysis has revealed the existence of 
a new instability, based on the cycle of vortical and 
acoustic perturbations in the subsonic part of the flow. 
The analytical study of this instability at high frequency 
(c/r s h «w« c/r son ) proves that it is independent of the 



postshock acceleration criterion established by Nakayama 
(1992, 1993). It is fed by the efficient production of vor- 
ticity perturbations when the shock is perturbed non ra- 
dially (Eq. ^2|), and by the vortical-acoustic coupling in 
the region of the sonic radius which enables the acoustic 
feedback. In this sense this non radial instability is generic 
and can be expected in more complex situations such as 
shocked flows with a weak angular momentum accreting 
into a black hole, even if the flow is decelerated immedi- 
ately after the shock. 

In the shocked Bondi flow the vortical-acoustic instability 
is faster than the radial instability if the shock is strong, 
by a factor oc log(A4i). 

The vortical-acoustic instability occurs for low degree per- 
turbations on a wide range of frequencies below the cut-off 
frequency ~ c/r son . A branch of unstable I = 1 eigen- 
modes corresponds to a vortical-acoustic cycle in which 
the acoustic feedback is produced in the pseudosound do- 
main (fshAsh < u < c / r sh)- The resulting growth rate 
is comparable to the growth rate of the vortical-acoustic 
cycle at high frequency. 

The role of the purely acoustic cycle was pointed out at 
high frequency in order to explain the large variations of 
the growth rate from one mode to another. More generally, 
the formalism developped in Appendix G concerning the 
simultaneous acoustic and vortical-acoustic cycles applies 
to any context where the efficiencies Q,1Z and timescales 
T Qi T n can be defined. 

- A strong I — 1 oscillatory instability was found in the 
pseudosound domain, at a frequency which is intermediate 
between the acoustic and advection frequencies (wf lax oc 

1/3 

A4 ± c/rsh). With comparable real and imaginary parts of 
the eigenfrequency, vorticity perturbations are advected 
over a very short distance during one growth time. On 
the basis of the contribution of a vortex to the Bernoulli 
perturbation sketched in Fig. |9[ this strong instability is 
a non radial consequence of post-shock acceleration. 

- As outlined in Sect. ||, the vortical-acoustic mecha- 
nism can be used as a tool in order to analyse the instabil- 
ity observed in numerical simulations of more complicated 
accretion flows. The specific application to BHL accretion 
or shocked flows with low angular momentum will be de- 
velopped in a forthcoming paper. 



Appendix A: Vorticity perturbations produced 
by the perturbed shock 

The non radial perturbation of velocity is deduced from 
the continuity of the velocity parallel to the shock, as in 
Landau & Lifshitz (1987, Chapter 90, p. 336): 



5v. 
6 u, 



= r d6 ' 
v\ — V2 <9A£ 
r sin 9 d(p 



(A.l) 
(A.2) 



14 



T. Foglizzo: vortical-acoustic cycle vs post-shock acceleration 



The perturbed vorticity in the flow is deduced from the 
Eqs. ( |A. l| ~ |A.2| ) and the Euler equation: 



we = (1-MiY 1- 



i-qcM 2 \ M\ dAv 



w v = -(1-M 2 2 f ( 1- 



wr s h / r sin 9 dtp 
ir]cM2\ M\ dAv 



06 



(A.3) 
(A.4) 



(A.5) 



The conserved quantity 5K s h defined by Eq. (i2|) is de- 
duced from Eqs. (A.4- A.5| ), resulting in Eq. (^0). 



Appendix B: Approximations of the acoustic 
perturbation 

B.0.1. High frequency limit: WKB approximation 



Appendix C: Formulation of the boundary value 
problem 



The differential system (g[j9|) is written as a single differ- 
ential equation of second order: 



dr 2 



d log fl- M' 

dr 



M 



2iuM 



M 



L 2 

rLf 
2 J 



SK. 



(l-M 2 )c 



df 
di- 



sh 



(1-M 2 )r 2 



(C.l) 



Follo wing the method used in F01, the general solution of 
Eq. (C.l) can be written using the solutions /o,/i of the 
homogeneous equation, where fo is the unique solution 
which is regular at the sonic radius. Let us normalize /q 
using the following definition of the WKB solutions f^ 
(Eq. B.l). The complex coefficient TV is defined by: 



If the wavelength of the perturbation is shorter than the 
lengthscale of the flow inhomogeneity (ojr/c ^ 1), a WKB 
approximation enables us to describe the propagation of 
acoustic waves in the direction of the flow (/q~) or against 
it(/o + ): 



ft 



M?c 2 



■ exp ilu 



1-M 2 ~'' 



(B.l) 



where we have chosen the normalization of f such that 
the lower bound of the integral is the turning point r t . 
The condition that r t < r s h defines a minimum frequency 
wf^ using Eq. (|ll|) and L 2 = 1(1 + 1): 



w />i 



w/(V S h) 



(B.2) 



B.0.2. Low frequency limit: uniform steady medium 
approximation 

Far from the accretor, the flow velocity decreases like cx 
1/r 2 and the density of the gas is uniform. If the shock is 
strong (r s h 3> r SQn ), the homogeneous equation associated 
with Eq. ((11) is approximated for r 3> r son and to <C 
cjirM) as follows: 



d 2 f 2df /to 2 L 2 . 
or z r dr \ c z r z 



(B.3) 



which is nothing more than the equation of acoustic waves 
in a uniform steady medium, in spherical coordinates. The 
solution / can therefore be approximated with a spherical 
Bcsscl function of the first kind ji, which is normalized 
here as in Eq. (B.l): 



Jo 



c \2ujrJ '+2 V c / 



3 



2 ^^son . 

e 4 c ji 

c \ c 



/ 7TC 

I 

(t) 



(B.4) 
(B.5) 



f = f++K*f . 



(C.2) 



The wronskien of the couple of solutions /q~, f Q is deduced 



from Eqs. (B.l) and (C.l) 



dr 



df 



h g r h - l 



2iuoMc: 



3 2i, 



J rt 1-M 2 



M 2 



(C.3) 



Let us normalize f± such that the wrons kien of the couple 
of solutions /o, /i is the same as in Eq. ( |C.3[ ). The general 
solution is then: 



SKsh -it» 

Q 

2iwc 3 



J r sh 



fo 



- fi 



i Mr 2 

fo 
Mr 2 



t _to f l+A^dr 



dr 



rr 1 + 

J rt 1- 



Mf dr 

Ta 1 m j r 



(C.4) 



A Frobenius analysis of /o,/i near the sonic points leads 
to: 



fo ~ /o(r s h) + 0{r - r Bon ), 
fl oc (r- r son )~M^. 



(C.5) 
(C.6) 



thus the integrals in Eq. (C.4) are converging when 
T r son if > c(dM/dr)(r son ). The regularity at the 
sonic radius therefore requires A\ = 0. A combination of 
Eq. (C.4) and its derivative at the shock radius leads to 
eliminate Aq as follows: 



dfo , df 
-^—J — E~Io 
dr dr 



M 



, .9 ^sh 



fo si i 



+ .M^ dr 
2" M 



Mr 2 



dr. 



(C.7) 



Eqs. ( p^[]I^ ) and ( p(i| ) provide the boundary conditions at 
the sho ck in order to replace f(r s h), df/dr{r s ^) and <5A" s h 
in Eq. flC.Tj) , resulting in Eq. (g|). 
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Appendix D: Calculation of Q a dv 

The analytical expression for Q a d v can be determined by 
writting the solution corresponding to zero acoustic flux 
F + at an outer boundary R, using the couple of homoge- 
neous solutions (Jq, f ~): 



f 



5K R 


e Jr * 


2iujc 3 








{/o 4 




B + f 






Jr Bon 






r R 


fo 





fo 



. r ii 

j„ i- 



Mr 2 



•i- ■ v ' 



dr* 



The regularity at the sonic radius requires 
B 1 = -TZ*B . 



(D.l) 



(D.2) 



The condition of absence of an incoming acoustic flux at 
the outer boundary is obtained by canceling the coefficient 
of /q 1 " at r = R, in the WKB limit of high frequency 
(w > c/R): 



Bo 



1 M 



Mr 



(D.3) 



Together with Eqs. (|D.2|) and 
boundary becomes: 



|),Eq. JD.ID at the outer 



f(R) 



e r t 

R 



1 /o 
/o ' Mr 2 " 



c Jr t 1 



+ .M^ dr 

T "dr. 



(D.4) 



Qadv(R), defined by Eq. (|2q), is deduced from the 
asy mptot ic behaviour of f$ in the WKB approximation 
(Eq. B.l). This calculation is formally similar to the cal- 
culation of Qk in F01, corrected for a phase shift. The 
efficiency Q a dv involved in the vortical-acoustic cycle is 
deduced from Eq. (|29|), with R — r^. In the strong shock 
limit, below the cut-off frequency (c/r son 3> u> ^> c/r s h), 
the acoustic efficiency Q a dv is approximated using the 
Spherical Bessel function ji (|B.5|) : 



is. 



adv 



A 



wr so „ 
c 



-i\x* 



ji(x)dx 



(D.5) 
(D.6) 



For A ^> 1, a function H(l) may be defined such that 

1 H(l) 



-i\x* 



ji(x)dx 



1 • 3. . . (2Z+ 1) A^ 



(D.7) 



The main contribution to the integral in Eq. (D.7) comes 
from the region x ~ A -1 / 3 . Eq. ( |30| ) is obtained from 
Eqs. and (O) 



Appendix E: Decomposition of the perturbation 
onto vorticity waves and acoustic 
waves in the WKB approximation 

The perturbations /, g immediately after the shock, de- 
fined by Eqs. (U|) and (|l^) are decomposed as follows: 



f(r ah ) = f_ + f + + f K , 
9{r sh) = .9- + .9+ +9k, 



(E.l) 
(E.2) 



where /k,9k correspond to the vorticity wave associated 
with the vorticity perturbation SK, and f±,g± correspond 
to the purely acoustic waves propagating in the direction 
of the flow (index +) or against the flow (index — ). An 
exact calculation can be made in the case of the reflexion 
of an acoustic wave Sp- with wavevector (k»,k±) on a 
plane shock in cartesian coordinates, in the absence of a 
gradient of M. The vorticity wave !k-,9k is advected at 
the velocity of the fluid: 

df K 



or v 



da 



K 



-9k- 



dr v 
Replacing these derivatives in Eqs. 



J), we obtain: 



A" 



9K 



M\c 2 SK 



v 2 L 21 



SK 



v 2 L 2 



9±=± 



M 2 c 2 
LE 



Eqs. ( |E.5[ ), flE.6| ) andj 
order to obtain Eq. (| 

M 2 c 2 Av 



f± = 



-(1 



) and 



(E.3) 
(E.4) 

(E.5) 
(E.6) 
(E.7) 

) are used with Eqs. in 



2 ,(fJ-TM 2 ) 2 



(E.. 



2fj, v y Al 1thM 2 

where fi in cartesian coordinates is also defined by 
Eq. (10), but replacing 1(1 + l)/r 2 by the wavenumber Ajj 



According to the definitions of Q s h and TZ s h in Eqs. (37) 
and (||), together with Eqs. <M) and (IeT" 

/ I ----- / ; :\/\ ■ > \ 

(E.9) 



li-M 2 \ fl + fiM 2 



[i + M 2 



Qsh 



-21(1 



1)1-^ 

'\M 2 



1 - iiM 2i 

1 (1 + fxM 2 )(l - Mj) 



(E.10) 



(V + M 2 ) 2 

These equations show the decrease of the vortical-acoustic 
coupling for weak shocks (Ai 2 ~ 1), and the existence of 
maximal efficiency at low frequency. Indeed, the maximum 

3 

Qsh ~ ^g-l(l + is reached for a frequency such that 
H ~ A4 2 /3, with K sh 1/4. 

The angle ip between the direction of propagation of the 
wave and the vector orthogonal to the shock surface is 
given by: 

UJ 
k I < 



tan-0 



M 2 



1 



(E.ll) 
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As remarked by Kontorovich (1958), the reflected sound 
wave propagates away from the shock with the same angle 
as the incident wave (?/>+ = "0-), as in a classical reflexion. 
The presence of a gradient of yV/f 2 in the Bondi flow pre- 
cludes the use of these formulae at low frequency. In the 
following calculation we assume ujr/c 3> 1 and .M2 <C 1, 
and keep only the first order terms in M 2 and c/cor. Thus 
[i ~ 1. Neglecting the coupling between the vorticity and 
acoustic waves in the vicinity of the shock, the vorticity 
wave fx, 9 k is s till advected at the velocity of the fluid. 
Eqs. flE.q ) and ( E.6 ) are now approximated at high fre- 
quency by 



9k 



M 2 2 c 2 5K 
5K 



(E.12) 
(E.13) 



Acoustic waves are described by Eqs. (|[]|) in the absence 
of vorticity perturbations, i.e. when SK = 0. Using the 



WKB approximation of Eq. (B.l), the radial derivative of 
f± is approximated by: 



df± 
dr 



g± is deduced from Eqs. (N) and (E.14) 



9± 



±1 



jyc \ f± 

2lut J M2C 2 



(E.14) 



(E.15) 



Eqs. flE.12| ), ( |E.13D and (|E.15| ) are used with Eqs. (|||^) 
in order to obtain Eq. (|20|) and 



M 2 c 2 Av / .r)c 

f± = {±l-M 2 -i— 

I v \ tor 



(E.16) 



Combining Eq. ( |20| ) with Eq. (|E.16[) , we obtain the ex- 
pressions for T^sh, Qsh in Eqs. (|40|) and (|64|). The pressure 
perturbations Sp± associated with th e acou stic waves f± 
are deduced from Eqs. (^0) and Eq. ( E.15| ): 



5p± = f± - v 2 g± 

p c 2 — V 2 

~ (1 T A4 2 )%. 



(E.17) 
(E.18) 



Appendix F: WKB approximation of the 
dispersion relation 

Using the integral expression of Q a dv (Eq. Eq. ( p2| ) 
can be approximated as follows for strong shocks (M2 <C 
1): 



77 - M2 



c 9/o 

2/(/+ 1) 
Ml 



+ 



iwr sh 



/o 



A4 2 ?7 



Qadv/o ■ 



(F.l) 



In the WKB limit ujr sh /c ^> 1, Eq. (F.l) can be simplified 
into 



/ + + ft*/-~-2Z(Z + l)Q. 



aclv 



_/o_ 

m| 



(F.2) 



Using the expressions for Q s h and 7?. s h (E qs. |40| and p^) , 
and the definition of Q and 1Z (Eqs. ^ and |58[),~hc global 
dispersion relation ( |49| ) is recovered. 

Appendix G: Analysis of the phase relation of the 
two cycles 

Let us introduce the new complex variable 

z=Qe luJTQ . (G.l) 

The resolution of the dispersion relation ([l9| ) is then equiv- 
alent to finding the complex number z satisfying: 



■ az 



1, 



(G.2) 



where the dimensionless parameter e = t-r/tq < 1 and 
I a I < 1 is defined by Eq. (p2]). The minimum and max- 
imum values z± of \z\ satisfy an equation similar to 
Eq. (E3) of FT00: 



z± =F \a\z± 



1. 



(G.3) 



The minimum and maximum effect of the acoustic cycle 
on t he gr owth rate uif 1 are directly related to z± through 



Eq. flG.lp : 

„ fs -L,„ g M. 

From Eq. (|G~3| ), 

z+ + z- = 2, 



(G.4) 



(G.5) 



we deduce that the maximum stabilizing effect of the 
acoustic cycle is to divide the efficiency Q by a factor 
2: 



1 < z+ < 2. 



(G.6) 



If the acoustic time is much shorter than the advection 
time (e <C 1), we deduce from Eq. (p.3|) the values of x±: 



z± 



1± lal 



(G.7) 



Thus the acoustic cycle may participate efficiently to the 
instability if e -c 1 and \a\ ~ 1. In the more general case 
where < e < 1, the following bounds on x are obtained 



from Eq. (|G_2|) 
l-la| < \z\ < 



1 



Q 



(G.8) 



The resolution of Eq. ( |G.2[) can be decomposed into a 
phase condition applied to the points of a continuous curve 
C, defined by: 



1| = \a\\z\ 



(G.9) 



T. Foglizzo: vortical-acoustic cycle vs post-shock acceleration 



17 



C is a closed curve containing the point (1, 0) in its interior, 
and (0, 0) in its exterior. It can be described in a univoque 
way by the angle tp = Arg(z — 1) 6 [— tt, tt]. Let us define 
the angle 9 = Arg(z) e [— tt, tt]. The solutions of Eq. (|G.2[ ) 
are recovered by applying to the solutions of Eq. ( p.9| ) the 
following phase condition: 



2kn = e( 



2k' tt), 



(G.10) 



where k, k' are two integers. Since (0, 0) is exterior to C, 
the range of values covered by 9 when ip covers [— tt, tt] 
is limited to \9\ < # max < tt. The discrete solutions of 



the phase equation (G.10) can be seen in a graphic way 
as the intersection of the periodic curve 9(tp) with the 
straight line 9 — <p/e. The number of solutions in each 
phase <p £ [0, 2tt] is therefore equal to 1/e on average. 
Comparing Eq. (G.l), Eq. ( p.S| ) with the definitions of 
<p, 9, one period of <p corresponds to a variation of lu t . of 
2tt/t-ji, and one period of 9 corresponds to 2tt/tq. The val- 
ues of |Q| and \TZ\ in Eqs. ( p4| ) and ( p5| ) can be determined 
from the measurement of uf^, and the average number 
n = tq/t-jz of eigenmodes per period Atu r = 2tt/tji, by 
eliminating |a|,at± from the set of equations (G.3), (G.4), 
and (|52 



Appendix H: Pseudosound approximation of the 
dispersion relation 

The integral involved in Eq. (|2^) can be approximated in 
the low frequency limit by introducing the complex vari- 
able z: 



dz _ . lu 1 + M 2 
dr ~ Mc 1 - M 2 ' 



(H.l) 



The contour of integration is deformed in the complex 
plane introducing a point on the real axis at +oo, and 
performing two integrations by parts: 



r s h \ 

J r son 




Mr 2 


r s hC , , 
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UT 


iuj 
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l + M* dr 
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fol-M 2 _ 



r 2 1 + M' 



e" z dz+ 



dz 



c 1-M 2 
iwr sh 1 + M 2 



fo + 



Mcr 2 h d ( f l-M 



iuj dr \ r 2 1 + M 2 



(H.2) 



The turning point in the Bondi flow would be at r t ~ 
Lc/ui. For r ^> r son , the homogeneous solution /o is ap- 
proximated by a Spherical Bessel function of the first kind 
jl(u>r/c). For ujr/c 1, it is approximated as follows: 



fo(r) ~ /o(r s h) 



r s h 



(H.3) 



The first integral on the right hand side of Eq. (H.2) is 
approximated using a Gamma function: 



r s hc 

ILU 



f l-M 2 { 
r 2 1 + M 2< 



z dz 



3Mc 
fo{r s h) 



(H.4) 



3Mz, 



+ 00 



sh 



/o(^sh) 



3MzJ 



-F 

V 3 



l + l 



(H.6) 



The second integral on the right hand side of Eq. (H.2) 
is negligible if r t 3> r s h 3> r 3on . In view of the particular 
case 1 = 1, the differential equation satisfied by /o is used 
to sum up terms of same order: 

Vr^ + l(l + l)fo = r 2 (l-M 2 ) d -^ + (^) 2 f 



dr 



\ c / 



2ivr\ df m „s 



Using Eqs. ( H.2 ), (H.6) and ( H.7 ), the dispersion relation 
( p2| ) is approximated in the pseudosound domain as fol- 
lows: 



2 d 2 f , 1(1 + 1) 



dr 2 



3z s h 



r— I fa dl ° g (— 
dr d log r \r 2 



3/^F 



/ + 4 



(H. 



The left hand side of Eq. (H.8) is dominated by the second 
derivative of /o if I > 2: 



l-l 1=2. 



''sh 



1 + 4: 

3 



(H.9) 



If I = 1, the second derivative of /o is approximated as 
follows: 



X 2 ^l , 
dx 2 

Thus Eq. 

Wz! h - 



x 
3' 



x 
~~5 



(H.10) 
(H.ll) 



) becomes: 

d log f M 
-I ° ' 

d log i 



4e Zsh r( §].(H.12) 
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